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Abstract 

We explored the underlying mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation (cell 
type switchings) from landscape and flux perspectives. Lineage reprogramming is a new regenerative method to convert a 
matured cell into another cell including direct transdifferentiation without undergoing a pluripotent cell state and indirect 
transdifferentiation with an initial dedifferentiation-reversion (reprogramming) to a pluripotent cell state. Each cell type is 
quantified by a distinct valley on the potential landscape with higher probability. We investigated three driving forces for 
cell fate decision making: stochastic fluctuations, gene regulation and induction, which can lead to cell type switchings. We 
showed that under the driving forces the direct transdifferentiation process proceeds from a differentiated cell valley to 
another differentiated cell valley through either a distinct stable intermediate state or a certain series of unstable 
indeterminate states. The dedifferentiation process proceeds through a pluripotent cell state. Barrier height and the 
corresponding escape time from the valley on the landscape can be used to quantify the stability and efficiency of cell type 
switchings. We also uncovered the mechanisms of the underlying processes by quantifying the dominant biological paths 
of cell type switchings on the potential landscape. The dynamics of cell type switchings are determined by both landscape 
gradient and flux. The flux can lead to the deviations of the dominant biological paths for cell type switchings from the 
naively expected landscape gradient path. As a result, the corresponding dominant paths of cell type switchings are 
irreversible. We also classified the mechanisms of cell fate development from our landscape theory: super-critical pitchfork 
bifurcation, sub-critical pitchfork bifurcation, sub-critical pitchfork with two saddle-node bifurcation, and saddle-node 
bifurcation. Our model showed good agreements with the experiments. It provides a general framework to explore the 
mechanisms of differentiation, dedifferentiation, reprogramming and transdifferentiation. 
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introduction 

A pluripotent undifTerentiated cell can differentiate into types of 
differentiated cells. Each cell type has a specific regulated gene 
expression. Cellular differentiation is determined by the underly- 
ing gene regulatory network during the process of development, 
which leads the primary cell into its ultimate fate-a particular 
phenotype. Induced pluripotent stem (iPS) cells provide the 
opportunity to obtain pluripotent stem cells which potentially 
have therapeutic uses [1,2]. Recentiy many studies have been 
reported that one type of cells can be converted to another type of 
functional cells directly [3-7] . This is a big step forward in the cell 
biology since there is no need to create iPS cells first for cell type 
switching, skipping many intermediate steps. This direct repro- 
gramming technology is called the lineage reprogramming. Thus 
an adult cell can be reprogrammed directly to new cells as lineage 
switching. The lineage switching through direct transdifferentia- 
tion without going through the iPS state might be applied to 



regenerative medicine with less risk of cancer. However, it is still 
challenging to quantify the mechanisms of the differentiation, 
dedifferentiation, reprogramming and transdifferentiation [3-11]. 

The concept of "epigenetic landscape" was first introduced by 
Waddington in 1940s [12] The quantifications of the Waddington 
potential landscape for the process of cell differentiation have been 
explored recently [13—17]. Different valleys represent different cell 
phenotypes (cell fates) on the cell development potential landscape 
[13-17]. Waddington visualized the undifferentiated state as the 
local maximum and differentiated states as the local minimum on 
the landscape [12]. In our landscape picture, the undifferentiated 
state and differentiated state are both local minima in certain 
regions of the landscape. Undifferentiated state has relatively low 
expressions of differentiation mark genes while differentiated state 
has at least one high expressions of differentiation mark genes. In 
addition, Waddington believed the differentiation is a downhill 
process driven by the funneled landscape gradient. In our picture, 
the differentiation can occur with several different mechanisms. 
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Figure 1. The scheme, phase diagram and intrinsic potential landscape of cell type switchings. A: The scheme of dedifferentiation 
(including reprogramming and differentiation) and transdifferentiation. B: A model for the gene circuit for cell development. C: The phase diagram for 
the gene circuit with 5' = 0.5,« = 4,A-i = fo = 1-0. D: The cell fate landscape i^q obtained from the Hamilton-Jacobi equation versus a and x, and the 
phase diagram was drawn on the intrinsic potential landscape with stable states represented by black solid lines and unstable states represented by 
black dash line. The red dash lines represent the dedifferentiation(reprogramming) and redifferentiation process while the yellow solid lines 
represents the transdifferentiation process. (a\ =02 = 1-0, ^1=^2 = 10, 5' = 0.5,« = 4,/i:i =/r2 = 1-0.) 
doi:1 0.1 371/journal.pone.01 0521 e.gOOl 



through funneled landscape, through stochastic fluctuations and 
the probability fluxes even when the landscape is not fijnneled 
towards the difierentiated states, and through induction. 

For development and difierentiation system, we represent a cell 
as a chemical system having given genomic makeup, with each 
and every possible phenotype as apotential "state" [18,19]. This is 
very much analogous to the notion of a polypeptide, as a chemical 
molecule, can have many difierent possible "conformational 
states", although each individual protein molecule has only a 
particular state at a given moment in time. This chemical 
definition of "the system" is important. Imagine that proteins are 
defined only through biological functions; then diflerent confor- 
mations of a polypeptide will be considered as "difiFerent 
molecules." Then the notion of spontaneous conformational 
change would not make sense. Indeed, there are still cell biologists 
who think different cells from the same person as different cells; 
rather than as a "same chemical system in different states" [18,19]. 
The process of the cell development can be viewed as the system 
moving from one valley (primary or stem cell phenotype) through 
bifurcation to another valley (differentiated cell phenotype) on the 
potential landscape. And the transdifferentiation process can be 
viewed as the system escaping from one stable difierentiated valley 



to another difierentiated valley through certain paths on the 
potential landscape shown in Figure 1(A). The differentiated cells 
[Sa) can switch to another lineage cell type [Sb] through an explicit 
pluripotent stable state (Sc). Indirect transdifferentiation mecha- 
nism which requires an initial dedifierentiation step Sa^Sc^Sb 
shown in Figure 1(A). It illustrates a difierentiated cell {Sa) 
reprogrammed back to a pluripotent state [Sc] with less 
difierentiated, and then can be re-differentiated to another type 
of differentiated cell {Sb) [3,5,6]. This is a possible strategy of 
pluripotent lineage reprogramming while the enhancement of 
efficiency is required. The underlying process is a transdifferentia- 
tion involving a stepwise dedifferentiation. In addition to indirect 
transdifferentiation, there is another lineage reprogramming 
approach: the direct transdifferentiation mechanism as 
Sa^So^Sb shown in Figure 1(A). Direct transdifferentiation is 
a mechanism of converting one type of differentiated cells to 
another type of differentiated cells without undergoing through a 
pluripotent state or progenitor cell type. The differentiated cells 
{Sa) down regulate their own cell-specific genes {X^) and activate 
the target cell-specific genes {X2), thus they can switch to another 
lineage cell type {Sb) through an explicit intermediate stable state 
{So) or a series of indeterminate states [3-5,8,9] . In our study, the 
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intermediate state is defined as an intermediate stable state with 
low or medium pluripotency and having very low expressions of 
the differentiation mark genes, while a series of indeterminate 
states are defined as a series of unstable states with low or medium 
pluripotency and very low expressions of differentiation mark 
genes in the course of lineage switching. Sridharan et al [20] 
showed that partially reprogrammed cells as an intermediate stage 
of the reprogramming process can switch to the completely 
reprogrammed iPS state. Thus the states of partially repro- 
grammed cells may exist along the paths from a differentiated state 
Sa or Sb to iPS state Sc- The research by Mikkelsen [2 1] showed 
that partially reprogrammed cells can be trapped at a common 
intermediate state. Thus the states of partially reprogrammed cells 
may exist along the paths from a differentiated state Sa to another 
differentiated state Sb through an intermediate So or indetermi- 
nate states. These intermediate state and indeterminate states may 
have certain expressions of stem cell marker genes and thus can be 
viewed as partially reprogrammed cells. This is support(;d by the 
observation that fibroblast cells specific genes are efficiently 
silenced and the embryonic reprogramming is not fuUy induced 
in partially reprogrammed cells [20]. We believe that different 
experimental and environmental conditions can lead to quite 
different results and change the topological structure of the 
potential landscape [20,21]. The partially reprogrammed cells 
may be trapped in certain regions in the gene expression space. 

In this study, we term direct transdifiFerentiation as transdifiFer- 
entiation and indirect transdifferentiation requiring an initial 
dedifferentiation or rc'programming step as dediffcrcntiation. The 
goal of rcgenerati\ e medicine can potentially be realized through 
the processes of differentiation, dedifferentiation, reprogramming 
and transdifiFerentiation [4] . Here we use cell type switchings short 
for the terms "difiFerentiation, dedifferentiation, reprogramming, 
and transdifiFerentiation". Recent advances have shown that there 
are three possible driving forces for cell type switchings: (1) 
Slofhaslic Fluctuations . Cells choose their pathways of difiFerenti- 
ation stochastically in the process of development without 
apparent regards to environment or history [22]. Some studies 
in cell development reveal that intrinsic stochasticity is an 
important mechanism for development [22]. The extrinsic 
fluctuations are also expected to play a role in cell development. 
Thus the fluctuations can be a driving force for the processes of 
cell type switchings. (2) Gene Regulation, cell type switchings can 
be achieved by the change of regulation strengths of their lineage 
specific genes in many studies [6,8,9,14,15]. (3) Induction. Lineage 
specific cells can be reprogrammed to a pluripotent state through 
over-expressions of some defined transcription factors [23,24]. 
Transfection of certain cell specific genes into the primary cells, 
and over-expressions of the target lineage specific genes as well as 
certain stem cell-associated genes can induce the processes of cell 
type switchings. 

Given the three driving forces for cell fate decision making, it is 
still challenging on how to quantify the processes of cell t}pe 
switchings on the landscape, and how to connect them to 
experiments. These processes of cell type switchings are controlled 
by their underlying gene regulatory network. The lineage-specific 
transcription factors play a critical role in the processes of cell type 
switchings. In this study, we explored a simple cell differentiation 
network module with autoregulation and mutual antagonism 
between transcription factors (lineage-specific genes) [15,17], 
which exists in many cell difiFerentiation processes, shown in 
Figure 1(B). The lineage-specific genes can strongly instruct the 
cellular lineage choice. The circuit is composed of a pair of self 
activating autoregulation and mutual inhibiting cross-antagonism 
cell-specific genes Xi and X2 [15,17]. In iPSC or ESC (embryonic 



stem cell), pluripotent genes are often highly expressed, and most 
lineage related genes are off. However, there are examples of gene 
regulatory circuits with the same architecture in our study which 
control binary decisions at branch points of cell differentiation in 
multi-potent cells. Such mutual antagonism gene circuit modules 
(where the self activation can also be indirect) in binary branch 
points of cell lineage commitment can often be found. A lot of 
studies have explored the primed multipotent common myeloid 
progenitor (CMP) can difFerentiate to either myeloid cell or 
erythroid cell in blood cell formation by mutual antagonism 
interaction of transcription factor gene GATAl and PU.l shown 
in Figure 2(A) [25,26]. GATAl and PUA are both self-activated. 
In the genetic regulation of the inner cell mass/trophectoderm 
lineage decision, Oct4 represses expression of Cdx2, and Cdx2 
represses expression of Oct4 to allow the segregation of inner cell 
mass and trophectoderm lineages [27,28]. Oc?4 and Cdx2 axe 
mutual inhibited and self-activated [27,28] shown in Figure 2(A). 
In the genetic regulation of the epiblast/primitive endoderm 
Uneage decision, antagonism between Nanog and Gata6 results in 
segregation of primitive endoderm and epiblast within the inner 
cell mass [27,29,30] shown in Figure 2(A). Nanog and Gata6 are 
also both self-activated [30] . These three circuits all can be viewed 
as Xi and X2 in our network. 

We win study this key network module to uncover the 
underlying functional mechanisms of cell type switchings. The 
phase diagram in Figure 1(C) suggests that the system can have 
five different phase regions, each of which has different underlying 
landscapes with different distribution of valleys. Furthermore, we 
show how stochastic fluctuation, gene regulation and induction 
induce the cell type switchings. The potential landscape and flux 
both direct the processes of cell type switchings. Probability flux 
provide a curling force breaking the detailed balance and lead the 
biological paths of cell type switchings to be deviated from the 
paths obtained by steepest descent gradient of the landscape. The 
forward and backward paths of cell type switchings are irrevers- 
ible, without passing through the saddle point. Furthermore, the 
flux can become the main driving force for cell type switching 
when the landscape is not biased towards the specific processes 
[16,31]. Barrier height and dynamic transition speed are used to 
quantify the global stability of the landscape topography. The 
stability here represents the ability for a cell to stay at a certain cell 
type state against certain fluctuations. In practice, the fluctuations 
in some cases maybe small but never zero. We uncover and classify 
four mechanisms of cell type switchings: super-critical pitchfork 
bifurcation, sub-critical pitchfork bifurcation, sub-critical pitchfork 
with two saddle-node bifurcation, and saddle-node bifurcation. 

Results and Discussions 

I. The model of cell fate network 

We start with gene circuit module for typical differentiation. 
The gene r('gulatory circuit for cell fate decision has two mutual 
repression and self-activation Uneage-specific transcription factors: 
Xi and X2 shown in Figure 1(B). It is more complete to consider 
three or more gene system. But the challenge is that a network 
with more genes requires more parameters to describe and 
therefore much bigger search space to explore exhaustively for 
uncovering the underlying mechanisms. Furthermore, with more 
genes, it is more diflRcult to visualize the results. The two gene 
system we considered is the simplest to exhaustively and effectively 
explore the underlying mechanism in parameter space [15-17,25]. 
We would like to use this model to explore the basic underlying 
mechanisms. The dynamics of this circuit is described by a set of 
two-variable ordinary differential equations below, with the rate of 
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Inhibition Macrophage B cell Mesendodermal Ectodermal 

Act i vat i on 



Figure 2. The gene circuits of mutual antagonism and self activation. A: The interaction of FUA and GATAl in determining myeloid cell or 
erythroid cell, Oct4 and Cdx2 in determining inner cell mass or trophectoderm, Nanog and GATA6 in determining epiblast or primitive endoderm. B: 
Scheme for the gene circuit of B cell to macrophage conversion. The dashed lines indicate uncertainty. C: Scheme for the gene circuit in determining 
mesendodermal and ectodermal. 
doi:1 0.1 371 /journal.pone.01 0521 6.g002 



expression change for these two genes: 



a\x" 



+ ■ 



dx2 aix^ biS" 
"df^ S"+x'' S" + x", 



-k\X\ =F\ 

-k2X2=F2, 



(1) 



where a'i and A-2 are the time-dependent expressions of the two 
cell-specific transcription factors Xi and X2 [15,17,25]. Parameter 
«! and 02 are the self activation strength of the transcription 
factors Xi and X2 respectively, bi and 62 are the strength of the 
mutual repression for transcription factors Xi and X2 respectively. 
ki and k2 are the first-order degradation rate for Xi and X2 
respectively [15,17,25]. S represents the threshold (inflection 
point) of the sigmoidal functions, i.e., the minimum concentration 
needed for appreciable changes, and n is the Hill coefficient which 
represents the cooperativity of the regulatory binding and 
determines the steepness of the sigmoidal function. For simplicity, 
we do not include studies of all the different parameters of 5 and n 
in the main text. We included the studies in the supporting 
information. We show the phase diagrams for varying these 
parameters in Figure SI in File SI. We can see varying these 
parameters can also lead to bi-stable states or tri-stable states and 
also the phase transitions. In the main text, the parameters for HiU 
function and degradation rate for Xi and X2 are specified as: 
5 = 0.5,n = 4, and A: = A:i =k2 = 1.0 [15-17,25]. In this section, we 
assume the symmetric situation 0 = 01=02 and b = b\=b2. 
Although the values of parameters can be different in organisms 
under diflerent circumstance, the mathematical model here 
describes a simple yet representative motif gene circuit, and these 
values (5' = 0.5,n = 4,A; = A:i = k2 = 1.0) are used in many previous 
studies [15-17,25]. 

1. The phase of cell fate network. To explore the dynamics 
under different conditions mimicking by different choice of 
parameters, we showed the phase diagram in Figure 1(C). If we 
can keep the mutual repression strength b frxed and the self 
activation o at various levels mimicking the actual developmental 
process where expression levels of transcription factor change 



[15](e.g. The expression level of transcription factor KlfA can be 
viewed as the effective self activation a at various levels mimicking 
the actual developmental process [32]. Because KlfA is not 
required for the maintenance of undifferentiated state of ES cells 
[32]. Furthermore, the expression level of KlfA decreases 
gradually after induced differentiation [32].), the cells are attracted 
to diflerent differentiated and undifferentiated states. There are 
five regions in the parameter phase space in Figure 1(C). Region I 
with lower self activation a and mutual repression b has only one 
stable state So with lower equal levels of the expressions of two 
lineage specific genes X\ and X2 shown in Figure 1(A). This is an 
intermediate state phase with lower lineage specific genes in the 
process of transdifferentiation [4] . Region II with higher mutual 
repression b and lower self activation a has two stable states shown 
in Figure 1(A): Sa which represents the differentiated state with 
higher expression of X\ and lower expression of X2, Sb which 
represents another differentiated state with lower expression oi X\ 
and higher expression of X2. Region III with lower mutual 
repression b and relative higher self activation o has three states: 
Sa,Sb and So- Region IV with higher mutual repression b and self 
activation o has three states: Sa,Sb and Sc which represents a 
pluripotent state with medium equal expressions of Xi and X2 in 
the process of dediflferentiation which can also be viewed as the 
process of reprogramming. Region V with lower mutual 
repression b and higher self activation o has all the four stable 
states: Sa,Sb,Sc and So- 

By changing the parameters of self activation o and mutual 
repression b, we can induce the initial differentiated cell to another 
diflerentiated cell in region II through the region III or region I by 
transdifferentiation (the yellow solid line), or through the region IV 
by dedifferentiation (the red dash line). In regions II, III, IV and V, 
there also exist tansdifferentiation within each. We will explore the 
dynamics of gene regulatory network for cell fate decision making 
process resulted from three driving force of stochastic fluctuations, 
gene regulation and induction through the instructive changes in 
details via the corresponding landscape topography for cell 
development. 

2. Super-critical and sub-critical pitchfork bifurcation 
versus saddle-node bifurcation in cell fate network. We 
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explored the bifurcation for cell fate decision network for different 
conditions. When mutual repression regulation parameters 
b = b\=b2 increase with small self activation regulation 
a = a\ =a2 = 0.2, the phase diagram has a super-critical pitchfork 
bifurcation which is a second order phase transition [33,34] shown 
in Figure 3(A). The solid lines represent stable frKcd points while 
the dash lines represent unstable fixed points. We can see a stable 
state So becomes an unstable state and splits into a pair of new 
stable states Sa and Sg at the critical point [33,35]. As the self 
activation regulation strength a increases, the phase diagram 
changes to a new form of sub-critical pitchfork with two saddle- 
node bifurcation which is a first order phase transition shown in 
Figure 3(B) as a = ai =«2 = 0.5. The initial state So is mono stable 
at lower mutual repression b, then a pair of new stable states Sa 
and Sb (two saddle-node bifurcations) emerge at somewhere far 
away from the initial state Sq as mutual repression b increases. 
After the critical point of sub-critical pitchfork, the center initial 
stable state Sq at the c(;iiter becomes unstable, only the two new 
stable states Sa and Sb are left in the phase space. Super-critical 
pitchfork bifurcation represents a type of "second-order transi- 
tion" in physics [36]. The difference between super-critical 
pitchfork bifurcation and sub-critical pitchfork bifurcation is that: 
super-critical pitchfork bifurcation represents one stable equilib- 
rium splits into two stable equilibrium and a unstable equilibrium 
while sub-critical pitchfork bifurcation represents two unstable 
equilibrium and a stable equilibrium merge into an unstable 
equilibrium. Thus super-critical pitchfork bifurcation differs from 
the sub-critical one in that two new stable equilibrium Sa and Sb, 
when they appear, already have a significant distance away from 
the middle stable equilibrium Sq- But the two stable frKcd points 
and the two unstable fixed points in sub-critical pitchfork with two 
saddle-node bifurcations are both symmetric in x\—X2—b three 
dimensional space, while they are not symmetric in X\—h two 
dimensional space shown in Figure 3(B). These two bifurcations 
shown in Figure 3(A) and (B) are similar to the picture described in 
Waddington's epigenetic landscape [12]. 

The phase diagrams shown in Figure 4(A), Figure 5(A) and 
Figure 6(A) are saddle-node bifurcations. A saddle-node bifurca- 
tion denotes a collision and disappearance of two equilibria rather 
than a pitchfork bifurcation [33,3.5]. The saddle-node bifurcation 
is a first order phase transition [33,34] . We can see that the initial 
valley Sa does not split into new valleys as the description of 
Waddingtons epigenetic landscape (a pitchfork bifurcation) [35]. 
New valleys Sb and Sc or So are born at somewhere far from the 
existing valley Sa in the state space. It is anther way of creating or 
eliminating the valleys from the potential landscape besides a 
pitchfork bifurcation [35] . The cell moves to the new valley Sc or 
So and Sb in sequence under fluctuations since its own valley 
disappears in another saddle-node bifurcation. We have already 
explored another form of bifurcation for cell fate network as self 
activation a decreasing with b = \.Q'm our previous study [15,17]. 
The phase diagram was drawn on the intrinsic potential landscape 
as the black lines in Figure 1(D) which is a sub-critical pitchfork 
[33,34] at the phase transition point (a = 0.774). 

We would like to explore these mentioned non-equilibrium 
phase transition under fluctuations and gene regulation. We might 
monitor the expressions of the differentiation marker genes in time 
and obtain the correlation functions. The singularity of the self- 
correlation function indicates the first order phase transition 
(saddle-node bifurcation) and the continuity of that shows the 
second order phase transition [33,34]. Thus we might distinguish 
these mechanisms of cell t)'pe switchings. We wUl explore these 
mechanisms of four bifurcations through our potential landscape 
theory in details below. 



3. Intrinsic potential landscape. We obtained the intrinsic 
potential landscape (j)^ (see the section of Methods) with Lyapunov 
properties to rjuantify the global stability by solving the zero 
fluctuation limit Hamilton-Jacobi equation and the associated 
intrinsic flux velocity in the zero noise limit [3 7] . The population 
potential landscape of cell development can be obtained through 
the exploration of the underlying probability dynamics, by solving 
the Fokker-Planck diffusion equation (see the section of Methods) 
[15]. The population potential landscape U is related to steady 
state probability distribution P,,, through —InP^s under fluctua- 
tions. The intrinsic potential landscape is quantified at the zero 
noise limit while the population potential landscape is quantified 
under finite fluctuations. Both show the global properties of the 
cell developmental process. Although intrinsic potential landscape 
gives less information (only at zero noise limit) about the network 
than population potential landscape, it can be used to quantify the 
global stabihty due to its nature of being a Lyapunov function [37]. 
We can Ulustrate two-dimensional potential landscape (the 
coordinates x\ and X2) to one dimension. One dimensional cross 
section coordinate X links Sa side minimum through Sc middle 
minimum to Sg side minimum, x represents the gene expression 
levels, .\:<0 shows gene X\ is dominant while .\:>0 shows X2 is 
dominant. If the s(;lf acti\'ation strength a dc-creases relatively 
slowly, relative to gene regulation in development, the potential 
landscape can be viewed as a succession of one dimensional 
potential slice. Figure 1(D) shows the intrinsic potential landscape 
for normal cell differentiation development process from plurip- 
otent state (Sc) to differentiated states [Sa and Sb) and the 
pluripotent reprogramming process from differentiated states {Sa 
and Sb) to pluripotent state (5c). We can see the intrinsic potential 
landscape 0o can be used to quantify the Waddington's picture 
and has almost the same shape with the population potential 
landscape [15]. 

The red dash lines and the yellow solid line shown in 
Figure 1(D) schematically described the lineage reprogramming 
process: dedifferentiation and transdifferentiation, respectively. 
The dedifferentiation process shows that differentiated state Sa 
follows a step backward to a pluripotent state Sc and then is 
induced to re-dififerentiate to another differentiated state Sb- While 
the transdiffc'rentiation process shows that differentiated state Sa 
converts directiy to another differentiated state Sb through certain 
intermediate stable state or not. Much work has been done on 
Uncage reprogramming and progress has been made in manipu- 
lating the key regulator gene to convert cell lineages [3-6,8,9]. The 
understanding of the underlying mechanism is still challenging. 
We will discuss the possible mechanisms of these lineage 
reprogramming process in detail using this simple gene regulatory 
circuit. 

We can see that when self activation a = a\=a2 is strong with 
higher mutual repression b = b\=b2=\-Q, the valley of the central 
pluripotent state Sc is much deeper and the system is attracted to 
this valley shown in Figure 1(D). As the strength of self activation a 
decreases, the valleys of side diflferentiated attractors Sa and Sb 
become deeper while the central pluripotent state Sc becomes 
weaker. When the strength of self activation a approaching to 
zero, the central state Sc becomes a ridge and therefore it is not 
stable while the side states Sa and Sb become stable. This result of 
intrinsic potential landscape with global Lyapunov property of 
global stability shows the similar mechanism with the result 
obtained from exploring the population potential landscape [15- 
17]. 

In order to quantify the stability of each state from the potential 
landscape topography, we can apply barrier height to measure the 
relative weights between different stable states. We showed barrier 
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Figure 3. The dynamics of super-critical and sub-critical bifurcations for cell type switchings. A: The phase diagram for changing the 
parameter b = hi =62 with a = a\ =02 = 0.2. B: The phase diagram for changing the parameter b = h\ = 62 with a = a\ =02 = 0.5. C: The quantified 
dedifferentiation and differentiation landscape and pathways for continuous changing parameter h with a = a\=a2 = 0.2. D: The quantified 
dedifferentiation and differentiation landscape and pathways for continuous changing parameter h with a = a\ =02 = 0.5. 
doi:1 0.1 371 /journal.pone.01 0521 6.g003 



height of intrinsic potential landscape versus the strength of self 
activation a in Figure 7 A. We set ^<l>osA='l'os~'l'oA ^'^'^ 
^'I'os c = "/"os ~ '^OC; where ijios is the value of the intrinsic potential 
landscape at the saddle point between state Sa and state Sc, (I>oa 
represents the minimum value of the intrinsic potential landscape 
at differentiated state Sa while i^oc represents the value of that at 
pluripotent state Sc- Barrier height /^(j>Qsc clecreases as a 
decreases, and state Sc vanished after the phase transition critical 
point Oci =0.774, where the system transits from three stable states 
{Sa,Sb,Sc) to two stable states {Sa,Sb)- It imphes that the 
attraction of state Sc becomes shallower. Barrier height of Ai^q^^ 
increases first, then decreases. It shows that the attraction of the 
dilferentiated state Sa {Sb) becomes deeper first, then becomes 
weaker after another critical point around ac2 = 0.4. So differen- 
tiated states Sa and Sb at 0^2 are more stable. Figure 7B shows the 
intrinsic potential barrier height Ai^gs^ has positive correlation 
with the population potential barrier height BasA= Bas — Bua 
under the diffusion coefficient D = 0.004 and D = 0.02, where Bas 
is the value of the population potential landscape at the saddle 
point between state Sa and state Sc, Bua represents the minimum 



value of the population potential landscape at the differentiated 
state Sa- The mean first passage time (MFPT) is useful to 
characterize the global stability if stochastic fluctuations are the 
dominant source of noise since it measures how the system can 
globally communicate from one state to another. The intrinsic 
barrier height Ai^qs^ ^nd the corresponding MFPT have the 
correlation of T ~ exp{ABas a) shown in Figure 7(C) with diffusion 
coeflTicient D = 0.02. 

A cell is a non-equilibrium open system with exchanges of 
energy and information from the outside environment. This leads 
to dissipation which is determined by both potential landscape and 
flux. The dissipation can give another global physical character- 
ization of the non-equilibrium system. Non-equilibrium system 
dissipates both energy and entropy in steady state, where the 
entropy production rate is equal to heat dissipation rate. The heat 
dissipation rate is formulated as HDR= J F-Jc/x [13,37-40], 
which increases first then decreases as self activation a decreases as 
shown in Figure 7(D). This indicates that larger area of the 
dominant probability flux leads to more heat dissipation because 
the system needs to consume more energy [37]. The system 
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Figure 4. The dynamics of transdifferentiation undergoing an intermediate state. A: The phase diagram for decreasing at induced the 
differentiated state Sa to the other differentiated state Sb through the intermediate state So- (ao = 1.0, 6 = 0.2) B: The barrier heights of the 
population landscape versus the parameter ai. C: The quantified transdifferentiation landscape and pathways for continuous changing parameter ai. 
doi:1 0.1 371 /journal.pone.01 0521 6.g004 




Figure 5. The dynamics of transdifferentiation undergoing a series of unstable states. A: The phase diagram for decreasing ai induced the 
differentiated state Sa to the other differentiated state Sg- (ao = 10, 6 = 0.3) B: The barrier heights of the population landscape versus the parameter 
fli. C: The quantified transdifferentiation landscape and pathways for continuous changing parameter ai. 
doi:1 0.1 371/journal.pone.01 0521 6.g005 
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Figure 6. The dynamics of dedifferentiation undergoing a pluripotent state. A: The phase diagram for decreasing a\ induced the 
differentiated state Sa to the other differentiated state Sb through the pluripotent state Sc- (oo =4.0, ft = 0.5) B: The barrier heights of the population 
landscape versus the parameter a\. C: The quantified dedifferentiation landscape and pathways for continuous changing parameter a\. 
doi:1 0.1 371 /journal.pone.01 0521 6.g006 



consumes more energy in the process of the development with 
three dominant states while the system consumes less at the 
beginning of cell development and at the end of cell development 
with less states. The heat dissipation rate provides a global 
characterization of cell development. It is intimately related to the 
robustness of the underlying network. 

II. The mechanisms of cell type switchings 

1. Stochastic Fluctuations. The cell type switchings at a 
given stage of development with different symmetric self 
activation a at fixed mutual repression h. The stochastic or 
inductive cell development can often be influenced by the external 
environment. We showed the paths of state transitions in cell 
development on the intrinsic potential landscapes for different self 
activation a with fixed mutual repression 6=1.0 due to stochastic 
fluctuations shown in Figure 8. We can see the green lines 
represent the reprogramming or dedifferentiation paths from 
differentiated state Sa or Sb to pluripotent state Sc. while the red 
lines represent the differentiation paths from pluripotent state Sc 
to differentiated state Sa or Sb shown in Figure 8(A)(B) when self 
activation a is relative stronger and the system has three stable 
states. Its worth pointing out that a green path from differentiated 
state Sa to pluripotent state Sc connected to a red path from 
pluripotent state Sc to another differentiated state Sb can provide 
a possible mechanism of the process of dedifferentiation first and 
then redifferentiation shown in Figure 8(A)(B). We also showed 
that both the green and the red lines represent the transdiflfer- 
entiation paths from one differentiated state to another differen- 
tiated state shown in Figure 8(C)(D) when self activation a is 
relative weaker and the system has only two stable states, just as a 



toggle switch. The intrinsic flux velocity ((Ju/ Pss)d^o) represented 
by purple arrows are perpendicular to the negative gradient of 
intrinsic potential ( — V^g) represented by the white arrows in 
Figure 8 (see the section of Methods). 

The cell type switchings processes at a given stage of 
development with symmetric changing mutual repression /; 
while fixing self activation a. We considered the potential 
landscape changing under fluctuations with varying mutual 
repression parameter b = b\=b2 at a given state with frxed self 
activation fl = 0.9. Figure 9(A) shows the phase diagram for 
changing mutual repression strength b. We can see that when 
mutual repression strength b decreases below 6 < 0.222, a new 
stable state So emerges. This is an intermediate stable state 
between differentiated states Sa and Sb- There are lower 
expressions of gene X\ and X2 in state So- Dashed lines represent 
the saddle point between stable states. As mutual repression 
6 < 0.222, the system has all four states Sa,Sb,So and Sc- The 
fluctuations in the system can enable stochastic switching among 
the stable states. Note that smaller mutual repression strength b 
here represents larger repression effect since the parameter b is in 
the numerator of an inhibition term with a positive sign. Smaller b, 
that is larger repression, leads the system towards intermediate 
state So, while larger b which represents smaller repression effect 
leads the system towards pluripotent state Sc- 

Any given cell may take a completely different route back to 
their pluripotent state in principle. Certain sequence of stages can 
emerge in the process of cell type switchings [4] . In experiments, if 
there are several pathways, one can collect the statistics and fmd 
out the relative probabilities of each path, giving the quantification 
of the path weights. In modeling, path integral weights are 
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Figure 7. The barrier height, escape time and dissipation rate for different self activation strength a with mutual repression 
strength 6= 1.0 under fluctuations. A: The intrinsic barrier height Aij!i versus a. B: The intrinsic barrier height Ac^j,^.^ versus the population barrier 
height Bosa in U for £) = 0.004 and £) = 0.02. C: The escape time Inrj from the valley 5'^ versus the intrinsic barrier height Aij>QSA- D: The dissipation 
rate versus the decreasing parameter a. 
doi:1 0.1 371 /journal.pone.01 0521 6.g007 



calculated by the action of the system analogous to the classical 
mechanical systems which determine the likeliliood of one path 
versus the other. We often used the dominant paths with the 
largest weights to represent the major pathways. We showed four 
dominant biological paths on the corresponding population 
landscape with different mutual repression strength b = 0.l5 (A), 
5 = 0.22 (B), b = 0.25 (C) in Figure 10. These processes are 
fluctuation or induction induced transition. The purple hnes 
represent the paths from state Sb to state Sa while the black hnes 
represent the paths from state Sa to state Sb [15,37]. We can see 
there are two dominant paths with the same color for transdiflfer- 
entiation from a certain differentiated state to another differen- 
tiated state in each sub figures, one path is through intermediate 
state So while the other path is through pluripotent state Sc- We 
also found the two different colored development paths between 
each two states follow quite different routes. It is irreversible 



between the forward dedifferentiation and the backward dediffer- 
entiation paths through the pluripotent state Sc, and between the 
two transdifferentiation paths through intermediate state So or 
without an explicit intermediate state. This illustrates the 
irreversibility of the developmental paths which can be verified 
from the ongoing and future dynamical experiments. 

The path weight represents the probability of each route for cell 
type switchings. It can be used to predict the probability of 
different routes for cell type switchings. The path probability can 
be obtained by the action ^(x) for cell development (See methods 
for details). We labeled Apc the action of the path through state 
Sc, and Apo as the action of the path through state So- 
Figure 9(C) showed the logarithm of dedifferentiation path 
probability through state Sc divided that of transdifferentiation 
through state So decreases as mutual repression strength b 
becomes weaker. This showed that the dedifferentiation path 
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Figure 8. The paths of cell type switchings with different self activation strength a. The paths of differentiation (A,6)/ dedifferentiation (A,6) 
and transdifferentiation (C,D) for different a in zero-limit fluctuations on the intrinsic potential ijfiQ. Purple arrows represent the intrinsic flux velocity 
((J,5//'.,,)^^q) while the white arrows represent the negative gradient of intrinsic potential ( — V^Jq))- 
doi:1 0.1 371 /journal.pone.01 0521 6.g008 



probability through state Sc decreases or the transdifferentiation 
path probability through state So increases as mutual repression 
strength becomes weaker. 

The purple arrows represent the direction of the probability flux 
J while the black arrows represent the direction of the negative 
gradient of population potential landscape —VU shown in 
Figure 10. We can see the flux is almost perpendicular to the 
negative gradient of the population potential landscape [13,37]. 
The dynamics of transdifferentiation and dedifferentiation pro- 
cesses are determined by both gradient landscape and probability 
flux. Probability flux provides a curling force breaking the detailed 
balance, and leads the system to stay at the non-equilibrium state. 
The gradient force attracts the system into stable valleys. The 
potential landscape and flux both direct the processes of cell type 
switching. Flux can lead a system to move on even a relatively flat 
landscape, e.g., the limit cycle attractor, thus "flux-directed 
differentiation" and "down-hUl-directed differentiation (Wadding- 
ton)" both can occur in cell development. "down-hUl-directed 



differentiation (Waddington)" leads to the exponential waiting of 
barrier crossing while "flux-directed differentiation" gives a much 
more precise timing. Flux also can lead the biological paths of cell 
type switchings to be deviated from the paths obtained by steepest 
descent gradient, and the corresponding paths of cell type 
switchings are irreversible. We would like to point out additional 
flux can emerge from epigenetics of slow (non-adiabatic) 
transcription and translation regulations [41] often encountered 
in eukaryotic cells. The flux generated by the slow time scales can 
lead to the new mechanism of differentiation and reprogramming 
[31,42]. The competition of barrier crossing and slow binding can 
lead to optimal speed of cell type switching. [31,42,43]. 

It is worth noting that even though state Sq disappears in 
Figure 10(C), there stiU exist transdifferentiation paths through a 
series of indeterminate states near (0,0) position. This provides the 
possible mechanism of two ways of lineage reprogramming. We 
labeled the saddle point between state Sa and state Sc as while 
the saddle point between state Sa and state So as sj. In 



PLOS ONE I www.plosone.org 



10 



August 2014 | Volume 9 | Issue 8 | e105216 



Mechanisms of Cell Type Switchings 




0.10 0.15 0.20 0.25 0.10 0.15 0.20 0.25 



b b 

Figure 9. The phase diagram, barrier height, probability of the dominant path and mean first passaging time for different mutual 
repression strength h. A: The phase diagram for changing mutual repression strength 6 = 61=62 with a = ai =07=0.9. B: The barrier heights 
versus the parameter h. C: The probability of the dominant path through the progenitor cell state Sc divided that of the path through the 
intermediate state So versus the inhibition strength b. D: The mean first passaging time through the two paths versus the inhibition strength b. 
doi:1 0.1 371 /journal.pone.01 0521 6.g009 




Figure 10. The flux on the population potential landscape. The flux on the population potential landscape with £) = 0.004. Purple arrows 
represent the flux (J„) while the black arrows represent the negative gradient of population potential landscape ( — Vf/)) for a = 0.9, 6 = 0.15 (A), 
6 = 0.22 (B), 6 = 0.25 (C). The black lines represent the pathways from state Sa to state Sb while the purple lines represent the pathways from state Sb 
to state Sa- 

doi:1 0.1 371 /journal.pone.01 0521 6.g01 0 
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Figure 9(B), we can see barrier height Bosio = Usi — Uo measur- 
ing the stability of intermediate state So increases and barrier 
height Bas2A = Usi — Ua measuring the degree of difficulty for 
transition from state Sa to state So decreases dramaticlly as 
mutual repression strength b decreases, where Ugi is the potential 
value at saddle point Sj and Uj is the minimum potential value at 
valley Sj. This implies that state So becomes more stable and 
robust as b decreases. 

We also can explore MFPT by T~exp(A(x)) [44]. Importantly, 
MFPT is also useful to characterize stability of the network for 
changing the regulations represented by the self activation a and 
mutual repression b under a small but fixed fluctuations (during 
the regulation changes or induction) mimicking the real environ- 
ments. Figure 9(D) showed MFPT along dedifferentiation and 
transdifferentiation paths versus mutual repression strength b. We 
can s(;c that tlu' transdifferentiation path through state So 
becomes more preferred than dedifferentiation path through state 
Sc, and MFPT becomes shorter for transdifferentiation path 
through state So as mutual repression strength decreases. In other 
words, transdifferentiation process is easier (harder) and the 
dedifferentiation process is harder (easier) when mutual repression 
is weaker (stronger). 

2. Gene Regulation, Decreasing self activation ai and 
increasing self activation 02 induce the transdifferentiation 
process fi-om state Sa to state Sb with lower mutual 
repression strength b. The instructive change of landscape 
via varying regulation strengths is another important mechanism 
in action for cell development. Down regulating the lineage 
specific gene for initial primary differentiated cell and up 
regulating the lineage specific gene for final target differentiated 
cell can induce transdifferentiation or dedifferentiation. We 
explored this mechanism below with changes in decreasing self 
activation ai for gene Xi and increasing self activation 02 for gene 

Self activation strength can be set for describing the time 
evolution of the self activation regulation parameters as: 
fli =ao exp( — .^1^) [25] which continuously decreases in time 
(down-regulates cell specific gene Xi for differentiated state Sa) 
and another self activation regulation strength 
^2 = flo[l ~6xp( — which continuously increases in time (up- 
regulates cell specific gene X2 for target differentiated state Sb) in 
cell developmental process due to the inffuences of the regulations 
of other genes. Xi and A2 are the rates for the decrease of self 
activations ai and 02- We assumed the same value of 
1 = ).\ =12 = lO^'' for simplicity for the latter calculations. At this 
value of X, self activation strength a\ and a2 decrease relatively 
slowly compared with regulation dynamics of gene X\ and X2. 
Thus the dynamics is a slow non-equilibrium relaxation process, 
flo is the scaled value of self activation a\ and a2 [25]. 

We explored the transdifferentiation mechanism below with 
decreasing self activation a\ and increasing self activation 02 with 
lower mutual repression strength b. Figure 4(A) shows the saddle- 
node bifurcation phase diagram for decreasing self activation 
strength a\ with lower b = 0.2 and smaller ao = 1.0. Figure 4(B) 
shows barrier height versus decreasing self activation ai with 
D = 0.005. We defined the saddle point between state Sa and state 
So as Si, and the saddle point between state Ss and state So as .S2. 
Barrier height is defined as: Ba^ij = Usi ~ Uj, where Usi is the 
potential \ alue of the i saddle point, and Uj is the minimum at 
valley Sj. Barrier height can quantify the degree of global robust 
and stability at a valley. We can see the cell stays at the monostable 
differentiated state Sa at the beginning of the transdifferentiation. 
As self activation ai decreases, an intermediate state So emerges. 
Valley Sa is much deeper than valley So due to barrier height 



BusiA of valley Sa being higher than that ot Ba^io- It means the 
differentiated state Sa is more preferred and more attractive than 
intermediate state So- The system is preferred to stay at state Sa 
with gene Xi being dominant. As self activation strength ai 
becomes weaker and self activation 02 becomes stronger, the valley 
of state Sa becomes shallower while the valley of state So becomes 
deeper. Then, the valley of state So is more attractive than that of 
state Sa since barrier height Ba^iA is lower than barrier height 
Busio, and gene Xi and X2 are both at lower expressions. After 
state Sa disappears, the cell is driven into intermediate state So- 
As s(4f acti\'ation strength fli decrease's furtlu'r, the other 
differentiated state Sg emerges, and barrier height Bas2B becomes 
higher than barrier height Basio- Finally, the cell is forced into 
state Sb- This process interprets the mechanism of transdiffer- 
entiation from state Sa to state Sb through an intermediate state 
So- 

The above results showed the dynamics at certain stage of 
transdifferentiation. We can sdso explore the continuous dynamics 
controlled by the set of equations below: 



dxi 
dt 



a\x; 



hiS" 



S"+x". S"+x"^ 



-k\x\ 



dx2 02X2 b2S" 
~dt ~ S'' + }^ S" + x^^ 

da\ 



— k2X2 



(2) 



dt 



= —ka\. 



where A =10 and a2 = ao — a\. The continuous time dynamics 
of down-regulating gene X\ and up-regulating gene X2 is shown in 

Figure 4(C) with o,, = 10,Z) = 0.005,/? = 0.2 using Eq.2. We 
obtained the transdifferentiation paths on the four dimensional 
potential landscape. The purple path is from state Sa to state Sb 
while the green path is the reverse transition both through 
intermediate state So- It implies that the system with small mutual 
repression strength b (large inhibition) prefers the transdifferentia- 
tion path through intermediate state So- Although transdiffer- 
entiation process does not seem to occur naturally, it has been 
observed in many experiments. For example, the exocrine cells in 
adult mice can transdifferentiate into jS-ceUs using defined factors 
for direct reprogramming without passing through a pluripotent 
state but through an unnatural intermediate state [4,8,9]. 

Figure 5(A) shows the phase diagram of saddle-node bifurcation 
under ao = l-0 and mutual repression /> = 0.3. Figure 5(B) shows 
barri(;r h(;ight versus self activation a\ with Z> = 0.005. We defined 
barrier height as Bust =Us— Ui, where Us is the potential value of 
saddle point s between state Sa and state Sb, and Ui is the 
potential value at state S,-. Figure 5(C) shows the paths and the 
landscape for continuous dynamics using Eq.2 with 
ao = l.O,D = 0.005,h = 0.3i. We can see the cell stays at differen- 
tiated state Sa with higher barrier height BOsa at first, then the 
landscape \'alley tilts the cell from state Sa to state Sb, barrier 
height BusB becomes higher than barrier BugA and the valley of 
state Sa eventually disappears. Finally, valley S^ becomes deeper. 
The weights of these two valleys exchange at the end of 
transdifferentiation process [35]. This process interprets the 
mechanism of transdifferentiation from stat(; S^ to state Sb 
directly without through a specific intermediate state but through a 
series of indeterminate states. This result can be used to explain 
the mechanism that the enforced expressions of CjEBP with 
endogenous PU- \ can reprogram B cell into macrophages [4,5]. B 
cell specific marker is CD\9 while the macrophage specific genes is 
Mac—\. The gene regulatory circuit is shown in Figure 2(B). B 
cell commitment factor Pax5 can up-regulate many B cell specific 
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genes (such as CZ) 19). The macrophage commitment factor 
C/EBP can up-regulate many macrophage cell specific genes 
(such as Mac— 1) and down-regulate B cell specific genes (such as 
Pax5) [4,5]. Transcription factor PU.l is needed in the proces.s of 
transdifferentiation. The gene PU.\ has the property of auto- 
activation. Mikkola's work indicated that C/EBP and Pax5 act in 
mutual antagonisms [5,45]. The dashed lines for the auto- 
activation indicate uncertainty in Figure 2(B). Thus we can reduce 
the gene regulatory circuit in to two markers of CZ) 19 and 
Mac — 1 similar as our mutual antagonistic and self activation Xi 
and X2 [4,5]. C/EBPs inhibit B cell commitment transcription 
factor (B cell-specific genes) which down-regulates B cell marker 
CD\9 {X\) in B cell, and co-activate macrophage specific genes 
which up-regulates its target marker Mac — 1 [X2) in macrophag- 
es. B cells pass through a series of indeterminate states with lower 
expressions of B cell-specific genes CD 1 9 [X\ ) and macrophage- 
specific genes Mac — 1 {X2), which does not seem to undergo an 
initial dediflferentiation [4,5]. 

Figure 11(A)(B) show the logarithms of MFPT versus barrier 
heights using the same parameters in Figure 4(B) and Figure 5(B) 
respectively. We can see the time spent from one state to another 
and barrier height have the relationship as: x ~ exp{Ba). It implies 
that the harder the system is out from one valley with higher 
barrier height, the longer the escape time is. 

We also explored the behavior for the system when regulation is 
not symmetrical, not only for the case when self-activation strength 
a\ is not equal to self-activation strength a2, but also for the case 
when self-activation strength 02 is not changing synchronously 
with self-activation strength a\. In Figure S2 in File SI, we showed 
the potential landscape of continuous dynamics with self-activation 
strength fl2 set as a constant (02 = 0.65) while the self-activation 
strength a\ continuously decreases. The other parameters are 
diffusion coefficient Z) = 0.005, mutual inhibition strength 
b\ =i>2 =0.2. We can see the cell may stay at difiFerentiated state 
Sa at first since the basin of differentiated state Sa is lower than 
differentiated state Sb when self-activation strength ai = 1.5, then 
the landscape basin tilts the cell from the differentiated state Sa to 



the intermediate state So, and the basin of state Sa eventually 
disappears. Finally, the basin Sb becomes deeper, and the system 
shifts from the intermediate state So to the differentiated state Sb- 
The green path is from state Sa to state Sb while the purple path is 
the reverse transition from state Sb to state Sa both through the 
intermediate state So- In Figure S3 in File SI, we showed the 
potential landscape of continuous dynamics with self-activation 
strength 02 set as a constant (02 = 0.1) while the self-activation 
strength a\ continuously decreases. The other parameters are 
diffusion coefficient £> = 0.005, mutual inhibition strength 
b\ =^'2 =0.2. We can see the cell may stay at differentiated state 
Sa at first when self-activation strength oi = 1 .0, then the cell shifts 
from the differentiated state Sa to the intermediate state So, and 
eventually the basin of state Sa disappears. The green path is from 
state >S^ to state Sq while the purple path is the reverse transition 
from state So to state Sa- Here, intermediate state So may 
represent the partially reprogrammed cells. 

Decreasing self activation a\ and increasing self 
activation 02 induce dedifiTerentiation process from state 
Sa to state Sb with higher mutual repression strength 
h. We assumed self activation a\ and fl2 at relatively higher 
average scaled values with flo=4.0 and relative larger mutual 
repression strength 6 = 0.5 to induce the initial cell undergoing 
through a balanced pluripotent state [24]. Figure 6(A) showed the 
saddle-node bifurcation phase diagram for decreasing self activa- 
tion strength a\ with flo = 4.0 at different time. Figure 6(B) showed 
the barrier height versus the parameter a\ with Z) = 0.005. We 
defined the barrier height as Bat = t/^,„„_v ^ Uj, where Ucmax is a 
constant relative maximum value of population potential land- 
scape and Ui is the minimum value of population potential at 
valley S,. We can see the system begins with a monostable 
differentiated state Sa with higher expression of ceU-specific gene 
Xi and lower expression of (cll-specific gene X2- As parameter ai 
decreases, a saddle node bifurcation emerges, giving rise to 
another differentiated state Sb with lower expression of cell- 
specific gene Xi and higher expression of ceU-specific gene X2- 
Initially, barrier height Boa of valley Sa is much higher than that 
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Figure 1 1 . Mean first passage time versus barrier height with different mutual repression strength b. A: The logarithm of the mean first 
passage time (MFPT) versus the barrier heights according to Figure 4(B). B: The logarithm of the mean first passage time (MFPT) versus the barrier 
heights according to Figure 5(B). 
doi:1 0.1 371/journal.pone.01 0521 6.g01 1 
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of valley Sb, thus valley Sa is much more stable than valley Sb- So 
the system prefers to stay at differentiated state Sa- As fli becomes 
weaker and the corresponding aj becomes stronger, two self 
activations 01,02 for two cell specific mutually exclusive genes 
Xi,X2 are over-expressing balanced (relative higher expression), 
another stable pluripotent state Sc with medium expressions of 
gene Xi and X2 emerges, and the potential landscape has three 
valleys. Valley Sa quantified by barrier height Boa is deeper than 
valley Sc and valley Sg quantified by barrier height Boc and Bag 
at the beginning of valley Sc emerging. As self activation Oi 
dccrc'ascs and 02 increases further, valley Sc and valley Sb 
become deeper while valley Sa becomes shallower. Barrier height 
Bac is higher than Boa and Bob at ai =02- Therefore, the system 
with differentiated state Sa shifts to under pluripotent state Sc as a 
process of dedifferentiation. A recent experimental studies [24] 
proposed a model for the coupled pluripotency module (self- 
activation of Oct4 and Sox2) and for the differentiation module 
with mutual antagonism between the MEs (mesendodermal) and 
ECTs (ectodermal) shown in Figure 2(C). MEs inhibit the 
activation between OctA and Sox2, then Oct4 can only activates 
gene MEs, and inhibits gene ECTs [24]. This process can be 
viewed as MEs have the effect of self activation. Thus, this module 
can be reduced to two mutucd antagonism gene MEs and ECTs 
with indirect self activation as our gene regulatory circuit of Xi 
and X2 shown in Figure 2(C). It implies that higher self activation 
strength fli and 02 being balanced can lead the differentiated cell 
back towards the pluripotent cell. As self activation Oi keeps on 
decreasing and 02 keeps on increasing, the valley of the other 
differentiated cell state Sb becomes deeper than that of pluripotent 
cell state due to barrier height Bob being higher than Bac and 
Ba^. Eventually, the valleys of Sa and Sc disappear at their 
saddle-node bifurcation [35]. Thus the cell leaves the pluripotent 
cell state Sc and is forced to enter into the other differentiated cell 
state Sb- The results showed the mechanism of dedifferentiation 
and redifferentiation. This mechanism can be used to explain 
many studies of cell dedifferentiation process during tissue 
regeneration both in vitro and in vivo [6]. For example, Pox5 is 
essential for initiating B cell commitment and is continuously 
required to maintain B cell lineage commitment [6,7,45]. Pax5 
deletion can convert committed B cells into hematopoietic 
progenitors with pluripotency [6,7,45]. It is partiy similar as the 
circuit in Figure 2(B) if we substitute Mac— I into other lineage 
specific genes. Pax5 deletion means down-regulating the B cell 
specific genes (such as BCs) as the effect of self activation ai. This 
gene regulation can lead B cells {Sa) to dedifferentiate to 
hematopoietic progenitors {Sc)- Then these cells can re-differen- 
tiate to T cell, macrophage or granulocyte {Sb) under appropriate 
culture conditions, such as the T-cell-deficient circumstance to 
reconstitute T cell development [6,7]. The appropriate culture 
conditions can be achieved by up-regulating the target cell genes 
as the effect of another self activation 02- 

The population potential landscape at different developmental 
stage of decreasing self activation parameter fl] after the relaxation 
process to a steady state among Xi and X2 is shown in Figure 6(C) 
using Eq.2. The green fine represents the dedifferentiated path 
from differentiated state Sa to another differentiated state Sb 
through pluripotent state Sc- The purple line represents the 
backwards dedifferentiated path from differentiated state Sb to 
another differentiated state Sa also through pluripotent state Sc- 
We can see the irreversible paths on the four dimensional 
population potential landscape due to non-zero flux. The 
dedifferentiated landscape and the paths can be quantitatively 
described for predictions. 



Decreasing mutual repression strength b induces 
differentiation and dedifferentiation process from state Sc 
to state Sa{Sb) with certain self activation a. Figure 3(A) 
shows the phase diagram of super-critical pitchfork bifurcation 
under self activation a = ai=a2 = 0.2 while changing mutual 
repression strength b- We can see the potential landscape of 
continuous dynamics shown in Figure 3(C) using Eq.2 with 
Z) = 0.005, self activation a = ai =02 = 0.2 and decreasing mutual 
repression b is similar to Waddington's epigenetic landscape 
[12,35]. A cell valley can form from an undiflferentiation state 
around So- So can be viewed as a stem ceil state* with lower 
expressions of differentiation gene markers while Sc can be viewed 
as the stem cell with medium expressions of the stem cell markers 
[17,46]. When decreasing mutual repression strength b, the initial 
valley splits into two other valleys and the initial valley becomes a 
ridge [35]. The cell will choose one valley as its fate. Figure 3(B) 
also shows the phase diagram of another form of sub-critical 
pitchfork with two saddle-node bifurcation under larger self 
activation strength a = ai =02 = 0.5 when increasing mutual 
repression b. The continuous potential landscape shown in 
Figure 3(D) using Eq.2 with Z) = 0.005, self activation 
a = ai =02 =0.5 and decreasing mutual repression b is also similar 
to Waddington's epigenetic landscape [12,35] except the sur- 
rounding of the critical point. Around the critical point, there 
coexist three stable states So, Sa and Sb- We also quantified the 
paths on the potential landscapes. The purple lines represent the 
differentiation paths from undiflferentiation state to differentiation 
state while the green lines represent the dedifferentiation or 
reprogramming paths. We can see the paths are irreversible even 
in the pitchfork bifurcation due to the existence of flux. This 
mechanism can describe the autonomous cell fate specification 
[47]. Stem cells must fulfill two tasks of self-renewal and 
generation of diflferentiated cells. In symmetric cell division, each 
stem cell can divide to generate either two daughter stem cells or 
two diflferentiated cells symmetrically while in asymmetric ceU 
division, each stem cell splits to one daughter stem cell and one 
daughter diflferentiated cell [48] . The pitchfork bifurcation in this 
study can represent an asymmetry event that a polarized mother 
stem cell splits into two daughter cell M or N with diflferent 
expressions oi Xi or X2- If daughter cell M with a vt.ry low value 
of Xi or X2, it might fall into diflferentiated state, while daughter 
cell N with relative higher expression of Xi or X2 still stays at the 
pluripotent state [35,48]. The asymmetric cell division usually 
occurs early in embryogenesis [35,47]. 

3. Induction of over expression. Cell fate is influenced by 
inductive stimulus from a group of surrounding cells [23,35]. 
Over-expressions of defined transcription factors can induce one 
ceU type to another ceU type which does not depend on gene 
regulations. This has been achieved in practice using over 
expression of stem cell marker transcription factors. In our 
previous gene circuit studies of cell fate decision making for stem 
cell differentiation and development [15,17], the two genes in the 
network are both differentiation markers. The idea is that the 
specific diflferentiation markers when imbalanced will give 
differentiation of one cell fate or the other (two side basins Sa 
and Sb)- A more balanced diflferentiation marker setup (between 
the two) will lead to iPS stem cell state (center basin Sc)- Our 
theoretical work [15,17] has predicted the possibility of the seasaw 
mechanism (balance or imbalance) of reprogramming. That is 
over-expressing both the concentrations of diflferentiation marker 
genes in a balanced way can induce and force diflferentiated cells 
into iPS .stem cells or pluripotent cells [15,17]. 

The two mutually exclusive differentiation markers MEs {Sa) 
and ECTs {Sb) shown in Figure 2(C) with balanced over- 
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expressions of key transcription linage specific factors can induce 
the lineage cell into pluripotent state (5c) instead of the stem cell 
markers OctA and Sox2 for pluripotency of reprogramming as a 
"seesaw model" [24]. Our theoretical work [15,17] has already 
predicted this possibility of expressing differentiation markers for 
reprogramming and the seasaw mechanism suggested in their 
work [24]. 

Significant efibrts have been made towards the experimental 
converting fibroblasts (Sa) to cardiomyocytes [Sb) by induction of 
over-expressing key genes. It is reported that direct transdiflfer- 
entiation can be achieved by over-expressing gene GataA, Thx5 
and Meflc from fibroblasts {Sa) to cardiomyocytes {Sb) [11]. 
GataA, Tbx5 and Meflc are core transcription factors during 
early heart development and can co-activate other cardiac gene 
expression [11]. So GataA, Tbx5 and Meflc can be viewed as 
cardiomyocyte cell specific gene CMs which have self activation. 
Over-expressing gene GataA, Tbx5 and Meflc {CMs) can 
transdififerentiate fibroblasts to cardiomyocytes not through a 
pluripotent state [11]. Another experiment showed that the 
indirect transdifferentiation can be achieved with an initial 
dedifFerentiation from fibroblasts {Sa) through pluripotent precur- 
sor-Cardiac progenitor {Sc) by over-expressing some stem cell 
markers Soxl,c — Myc,OctA and KlfA, and then be induced to 
cardiomyocytes {Sb) [10]. 

Conclusions 

In this study, we applied our potential and flux framework to 
explore the mechanisms of cell developmental processes of 
dilferentiation, dedifFerentiation, reprogramming and transdifier- 
entiation. The potential landscape of two gene regulatory circuit 
shows that the system has four stable valleys at specific regulation 
regions, two differentiated state Sa and Sb, one pluripotent state 
Sc, and an intermediate state So- Our work provides a 
quantitative basis for explaining the mechanisms of the transition 
among the four states. Barrier height based on the population 
potential landscape or the intrinsic potential landscape can 
quantify the stabiUty of the attractors and the eSiciency of 
switching among the attractors. We can acquire the dynamical 
transition rate of the system from one valley of attraction to 
another by MFPT for escape and the dominant paths for 
dedifFerentiation and transdilferentiation via the path integral 
method. We can see the paths of cell type switchings are 
irreversible due to non-zero probability flux. 

In this study, we have discussed three driving forces: stochastic 
fluctuations, gene regulation and induction, which can lead to cell 
type switchings. The cell type switching driven by stochastic 
fluctuations is a spontaneous transition, gene regulation is much 
like a non-autonomous varying of time-dependent landscape, and 
induction is a condition of initial value re-setting process with no 
apparent paths. The fluctuations maybe small in some cases but 
never zero. When exploring the stochasticity, we used fixed set of 
the values of self activation and mutual repression regulation 
parameters a and b. We not only discussed the possibiUty of cell 
type switching through stochastic dynamics but also other two 
mechanisms including the induction and regulation changes. We 
also explored the different dynamics with dilferent sets of the 
parameter a and b. For gene regulation, we varied the parameters 
a and h for regulating the cell type switchings. For induction, we 
did not change any parameters. Instead, we just gave the cell an 
initial set (condition) with over-expression of its lineage specific 
gene. We quantitatively investigated the mechanism of cell type 
switching through the induction without the change of the 
underlying landscape and through the changes in regulations 



leading to the changes of the underlying landscape topography. 
Furthermore, these two types of cell type switchings driven by gene 
regulation and induction are not spontaneous transitions only due 
to fluctuations, but a controlled process under either the changes 
in regulations with regulation-dependent potential landscape or 
the induction with fixed potential landscape. 

We found that the topography of the global potential landscapes 
is strongly correlated to the self activation strength and the mutual 
repression strength of the transcription factors. Dedifferentiation 
can be induced by the core regulators of pluripotent genes using in 
iPS or the synergistic effect of lineage specifiers in specification of 
difierentiated cells. We can adjust two self activation strength ai 
and 02 to be relative larger to force the differentiated cell to a 
pluripotent cell with higher inhibition strength b, and then re- 
diflerentiate the pluripotent cell to our target differentiated cell 
type [24]. This process can be viewed as an initial epigenetic 
activation phase representing the redifferentiation after a temporal 
overexpression of pluripotent reprogramming factors to a plurip- 
otent state [4,24,49]. Somatic cells can be transdifierentiated by 
temporal over-expressions of pluripotent reprogramming tran- 
scription factors. TransdiflFcrcntiation can be induced by down 
regulating the lineage specific marker gene {X\) of the original 
differentiated cell (decreasing self activation fli) while activating 
another Uncage specific marker gene {X2) of the final differentiated 
cell (increasing self activation ai) at relative lower inhibition 
strength b, through an intermediate state or a series of 
indeterminate states. This process can be viewed as lineage- 
instructive transcription representing the induction of lineage 
specific gene for the target differentiated cells [4,49]. This gives us 
a new understanding that the topography of underlying potc-ntial 
landscapes in cell development dynamics determines the feasibility 
and efficiency of cell type switchings. 

We also classified the mechanisms of pitchfork bifurcations 
depicted Waddington's epigenetic development landscape includ- 
ing super-critical pitchfork bifurcation, sub-critical pitchfork 
bifurcation, sub-critical pitchfork with two saddle-node bifurca- 
tion, and saddle-node bifurcation depicted the transdifferentiation 
landscape [35]. We uncovered a pitchfork bifurcation of 
Waddington's epigenetic landscape and the irreversible paths 
(caused by the non-equilibrium flux) between diEF(;r(;ritiation and 
reprogramming. We also uncovered the saddle-node bifurcation 
landscape. Saddle-node bifurcation can give the explanation of 
possible mechanisms of dedifFerentiation and transdifferentiation 
processes and can further explain the irreversibility of the paths for 
differentiation, dedifFerentiation, reprogramming and transdiflFer- 
entiation processes as hysteresis loop even without the presence of 
the non-equilibrium flux. We noticed a special kind of sub-critical 
pitchfork with two saddle-node bifurcations also shares the certain 
features with saddle-node bifurcation (hysteresis loop) and certain 
features of pitchfork bifurcation (Waddington's landscape). 

Importantiy, we uncovered some novel mechanisms as a 
starting point to decipher the mysterious code of the cell type 
switchings. Our theory can be used to guide the designs of the 
differentiation, dedifferentiation, reprogramming and transdiffer- 
entiation processes. 

Methods 

Quantifying non-equilibrium potential landscape, flux, 
non-equilibrium thermodynamics and the paths 

Fluctuations exist widely in biological systems [13,50-55]. The 
dynamics in noisy fluctuating environments can be formulated as: 
x = F(x) + f. F(x) is the deterministic force, where x is the vector 
representing difierent concentrations in state space. { is Gaussian 
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noise term and its autocorrelation function is 
<f(x,f)^(x,0)> =2D(x)(5(0, where D(x) is difiFusion coefficient 
matrix. We set D(x) = DG(x), where D is the diffusion coefficient 
representing the level of noise strength while G is the scaled 
diffusion matrix described the anisotropy phenomenon. We can 
explore the corresponding Fokker-Planck diffusion equation 
[56,57] for probability distribution P(x,0: dP/dt = 
- V (F/'-DG(x)VP). In this study, we set G(x) as a unit matrix 
for simplicity. The probability flux J is defined as: J = FP—DWP. 
In steady state, the force decomposition is shown as: 
F=-Z)V[/ + J,,/P,, [13,50,51]. 

We obtained the Lyapuntn- function (/iq as the intrinsic potential 
from the zero fluctuation limit HamUton-Jacobi equation(HJE) 
[37,58]: F-V^o + Wo'Ci'V^o =0 by a numerical method - level set 
method using the Mitchell's level-set toolbox [59]. The force 
decomposition in zero fluctuation limit is shown as: 
F= — G-V^o + (J.s.s/^.!.!)lz)^o- r"rom the Hamilton-Jacobian equa- 
tion above, we can obtain (Js.s/P.!.j)li)^o'V0o = 0 [37,51,60,61]. 
We also can obtain the mean first passage time from the following 
equation [56]: F-Vt-|-Z)V^t= — 1. represents the mean first 
passage time from state / to state / 

The path integral approach we used is shown as below. The 
path probabilit)' starts from initial state X, at t = 0, and end 
at the final state of Xf at time t. The path integral formula is 

shown as [15,44]: P(x/,?|x,-,0) = Dx exp[- dt{^V-¥{x) + 
^ (dx/dt - F(x))--^-(dx/dt-F(x)))] -- 



Dx exp[—A{x)] = 



4 

J Dx exp[— J L(x(t))dt], where IAx{t)) is the Lagrangian and 
^(x) is the action for each path [15,44]. The path integral over Dx 

represents the sum over all possible paths connecting x, at time 
t = Q to Xf at time t. The exponent factor gives the weight of each 



specific trajectory path. The probability from initial state to the 
final state is equal to the sum of all possible paths with different 
weights. Every dynamical path doesn't contribute to the same 
weight and each path is exponentially weighted. Therefore, the 
path integrals can be approximated with a set of dominant paths 
while the other subleading path weights can be neglected for their 
relative small values. We can find the dominant paths with the 
optimal weights through minimization of the action ^(x) or 
Lagrangian L{x{t)). Thus, we can identify the optimal paths which 
give more contribution to the weight as biological paths or cell 
type switching pathways in our study. 
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